## pipeline script for MethyPattern align in to genome feature
## author: Yaping Liu  lyping1986@gmail.com
my $mem="10g";
#my $data_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/MethyPatternFeatureWalker/GBM_hg19/MethyPatternResult/NDR/";
#my $data_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/MethyPatternFeatureWalker/GBM_hg19/MethyPatternResult/NDR_segment_80bp_criteria/";
my $data_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/MethyPatternFeatureWalker/GBM_hg19/MethyPatternResult/NDR_segment_distance/";
#my $feature_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/MethyPatternFeatureWalker/GBM_hg19/location_data/";
my $feature_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/MethyPatternFeatureWalker/GBM_hg19/location_data/NDR/";
my $plot_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/MethyPatternFeatureWalker/GBM_hg19/plot/NDR/";
#my $data_dir = "/export/uec-gs1/laird/shared/research/yaping/code/GNOME_seq/MethyPatternFeatureWalker/GBM_hg19/MethyPatternResult/promoter/";
#my $feature_dir = "/export/uec-gs1/laird/shared/research/yaping/code/GNOME_seq/MethyPatternFeatureWalker/GBM_hg19/location_data/";
#my $plot_dir = "/export/uec-gs1/laird/shared/research/yaping/code/GNOME_seq/MethyPatternFeatureWalker/GBM_hg19/plot/promoter/";
my $genome_assembly = "hg19";
my $bissnp_dir = "/home/uec-00/yapingli/software/BisSNP/BisSNP-0.71_forNOMeseq_v6.jar";
my $R_plot_script = "/home/uec-00/yapingli/code/mytools/R/MethyPatternFeaturePlotSinglePlot.R";

#my $bed_file = "/export/uec-gs1/laird/users/yaping/code/tmp/AdamB/data/GM_WGBS_CpGsites_10xCov.sorted.bed";
#my $bed_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/NDR/HMM/GBM_hg19/NDR_segment_minCT3/";
#my $bed_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/NDR/HMM/GBM_hg19/NDR_segment_80bp_criteria/";
my $bed_dir = "/export/uec-gs1/laird/users/yaping/code/NOMeseq/NDR/HMM/GBM_hg19/NDR_segment_4_category/";
my $reference_seq;
my $dbsnp;
if($genome_assembly eq "hg19"){
	$reference_seq = "/home/uec-00/yapingli/data/genome/hg19_rCRSchrm.fa";
}
else{
	$reference_seq = "/home/uec-00/yapingli/data/genome/Homo_sapiens_assembly18_without_random.sort.fasta";
}
#my $sample="NS-157";
#my $sample="GM_WGBS";

my $data_matrix_scale = 2000;
my $confidance_level = 20;
my $feature_name = 1;

my $plot_step = 20;
my $plot_scale = 1000;
my $plot_axis_step = 250;
my $smooth = 1; ## 0: no smooth; 1: smooth

my $count = 1;
my $use_age = "USAGE: perl MethyPatternFeaturePipelineByBed.pl data_dir feature_dir plot_dir";
#if($#ARGV <3){
#	print "$use_age\n";
#	exit(1);
#}
opendir(DH2,$bed_dir);
foreach my $bed_file(readdir(DH2)){
	next if $bed_file eq "." or $bed_file eq "..";
	next if $bed_file =~ /\.idx/;
	#next if $bed_file !~ /NS157_NDR_Cluster1.6column.sort.bed$/;
	opendir(DH,$feature_dir);
	foreach my $feature_file(readdir(DH)){
	if($feature_file =~ /(\S+)\.minusUpstream0\.plusDownstream0(\S+sort.bed)$/){
	#if($feature_file =~ /(NS157_NDR_Cluster3.6column.sort.bed)\.minusUpstream0\.plusDownstream0(\S+sort.bed)$/){
		my $prefix = $1;
		my $suffix = $2;
		my $region_file = $prefix."\.minusUpstream3000\.plusDownstream3000".$suffix;
		my $header = "#PBS -q lairdprio\n";
			$header .= "#PBS -l walltime=168:00:00,mem=$mem,nodes=1:ppn=8\n";
			#$header .= "#PBS -l walltime=168:00:00,mem=$mem\n";
			#$header .= "#PBS -W depend=afterany:1996998.hpc-pbs\n";
			$header .= "#PBS -N $feature_file-$bed_file\n";
			$header .= "#PBS -j oe\n";
			$header .= "cd \$PBS_O_WORKDIR\n";
			my $java_cmd = "java -Xmx$mem -jar $bissnp_dir -T MethyPatternFeatureByBed ";
			$java_cmd .= "-R $reference_seq -values $bed_dir"."$bed_file -orientated ";
			$java_cmd .= "-feature $feature_dir"."$feature_file -distance $data_matrix_scale -minCTdepth 1 -alignmentType Center ";
			$java_cmd = $java_cmd."-L $feature_dir"."$region_file ";	
			my $gch_file = $data_dir."$prefix.$genome_assembly.$bed_file.NDR_density.$data_matrix_scale.txt";
			$java_cmd = $java_cmd."-dataPoint $gch_file ";
			
			$header .= "$java_cmd\n";
			$prefix .= ".$bed_file";
			my $R_cmd_plot = "/home/uec-00/shared/production/software/R-2.14.1/R-2.14.1/bin/R --no-restore --no-save --args wd=$plot_dir prefix=$prefix gchfn=$gch_file step=$plot_step scale=$plot_scale axistep=$plot_axis_step smooth=$smooth < $R_plot_script \n";
			$header .= "$R_cmd_plot ";			
			my $outfileName = "submitPBS" . $count;
			open(OUT, ">$outfileName") or die "can not open file:$!";
			print OUT $header;
			close OUT;
			system("qsub $outfileName");
			$count++;
			print "$feature_file is submitted\n";
			unlink $outfileName;
			#print "$header\n";
	}
}
closedir(DH);
}

